import numpy as np

def F4_MPC_Matrices_PM(A, B, Q, R, S, N_P):
    n = A.shape[0]
    p = B.shape[1]

    Phi = np.zeros((N_P * n, n))
    Gamma = np.zeros((N_P * n, N_P * p))
    Delta = np.zeros((N_P * n, n))

    for i in range(N_P):
        Phi[i * n:(i + 1) * n, :] = np.linalg.matrix_power(A, i + 1)

    for i in range(N_P):
        for j in range(i + 1):
            Gamma[i * n:(i + 1) * n, j * p:(j + 1) * p] = np.linalg.matrix_power(A, i - j).dot(B)

    # 填充每一行
    for i in range(N_P):
        if i == 0:
            # 第一行是单位矩阵
            Delta[i*n:(i+1)*n, :] = np.eye(n)
        else:
            # 其余行是 A 的 i 次方
            Delta[i*n:(i+1)*n, :] = np.linalg.matrix_power(A, i)

    Omega = np.kron(np.eye(N_P), Q)
    Omega[-n:, -n:] = S

    Psi = np.kron(np.eye(N_P), R)

    M = Gamma.T.dot(Omega).dot(Phi)
    H = Psi + Gamma.T.dot(Omega).dot(Gamma)
    F = Gamma.T.dot(Omega).dot(Delta)
    return Phi, Gamma, Omega, Psi, M, H, F
